home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / LinAlg 3.1 / LinAlg / vmatrix1.cc < prev    next >
Text File  |  1995-12-20  |  10KB  |  329 lines
  1. // This may look like C code, but it is really -*- C++ -*-
  2. /*
  3.  ************************************************************************
  4.  *
  5.  *        Verify Advanced Operations on Matrices
  6.  *        (multiplication, inverse, determinant evaluation)
  7.  *
  8.  * $Id: vmatrix1.cc,v 3.1 1995/02/03 15:32:08 oleg Exp oleg $
  9.  *
  10.  ************************************************************************
  11.  */
  12.  
  13. #include "LinAlg.h"
  14. #include <math.h>
  15. #include "builtin.h"
  16. #include <iostream.h>
  17. #include <float.h>
  18.  
  19. /*
  20.  *------------------------------------------------------------------------
  21.  *            Verify the determinant evaluation
  22.  */
  23.  
  24. static void test_determinant(const int msize)
  25. {
  26.   cout << "\n---> Verify determinant evaluation\n"
  27.           "for a square matrix of size " << msize << "\n";
  28.  
  29.   Matrix m(msize,msize);
  30.  
  31.   cout << "\nCheck to see that the determinant of the unit matrix is one";
  32.   m.unit_matrix();
  33.   cout << "\n    determinant is " << m.determinant();
  34.   assert( m.determinant() == 1 );
  35.  
  36.   const double pattern = 2.5;
  37.   cout << "\nCheck the determinant for the matrix with " << pattern <<
  38.           "\n    at the diagonal";
  39.   register int i,j;
  40.   for(i=m.q_row_lwb(); i<=m.q_row_upb(); i++)
  41.     for(j=m.q_col_lwb(); j<=m.q_col_upb(); j++)
  42.       m(i,j) = ( i==j ? pattern : 0 );
  43.   cout << "\n    determinant is " << m.determinant();
  44.   assert( m.determinant() == pow(pattern,(long)m.q_nrows()) );
  45.  
  46.   cout << "\nCheck the determinant of the transposed matrix";
  47.   m.unit_matrix();
  48.   m(1,2) = pattern;
  49.   Matrix m_tran(Matrix::Transposed,m);
  50.   assert( !(m == m_tran) );
  51.   assert( m.determinant() == m_tran.determinant() );
  52.  
  53.   cout << "\nCheck the determinant for the matrix with " << pattern <<
  54.           "\n    at the anti-diagonal";
  55.   for(i=m.q_row_lwb(); i<=m.q_row_upb(); i++)
  56.     for(j=m.q_col_lwb(); j<=m.q_col_upb(); j++)
  57.       m(i,j) = ( i==(m.q_col_upb()+m.q_col_lwb()-j) ? pattern : 0 );
  58.   assert( m.determinant() == pow(pattern,(long)m.q_nrows()) * 
  59.      ( m.q_nrows()*(m.q_nrows()-1)/2 & 1 ? -1 : 1 ) );
  60.  
  61.   cout << "\nCheck the determinant for the singular matrix"
  62.           "\n    defined as above with zero first row";
  63.   m.clear();
  64.   for(i=m.q_row_lwb()+1; i<=m.q_row_upb(); i++)
  65.     for(j=m.q_col_lwb(); j<=m.q_col_upb(); j++)
  66.       m(i,j) = ( i==(m.q_col_upb()+m.q_col_lwb()-j) ? pattern : 0 );
  67.   cout << "\n    determinant is " << m.determinant();
  68.   assert( m.determinant() == 0 );
  69.  
  70.   cout << "\nCheck out the determinant of the Hilbert matrix";
  71.   Matrix H(3,3);
  72.   H.hilbert_matrix();
  73.   cout << "\n    3x3 Hilbert matrix: exact determinant 1/2160 ";
  74.   cout << "\n                              computed    1/"<< 1/H.determinant();
  75.  
  76.   H.resize_to(4,4);
  77.   H.hilbert_matrix();
  78.   cout << "\n    4x4 Hilbert matrix: exact determinant 1/6048000 ";
  79.   cout << "\n                              computed    1/"<< 1/H.determinant();
  80.  
  81.   H.resize_to(5,5);
  82.   H.hilbert_matrix();
  83.   cout << "\n    5x5 Hilbert matrix: exact determinant 3.749295e-12";
  84.   cout << "\n                              computed    "<< H.determinant();
  85.  
  86.   H.resize_to(7,7);
  87.   H.hilbert_matrix();
  88.   cout << "\n    7x7 Hilbert matrix: exact determinant 4.8358e-25";
  89.   cout << "\n                              computed    "<< H.determinant();
  90.  
  91.   H.resize_to(9,9);
  92.   H.hilbert_matrix();
  93.   cout << "\n    9x9 Hilbert matrix: exact determinant 9.72023e-43";
  94.   cout << "\n                              computed    "<< H.determinant();
  95.  
  96.   H.resize_to(10,10);
  97.   H.hilbert_matrix();
  98.   cout << "\n    10x10 Hilbert matrix: exact determinant 2.16418e-53";
  99.   cout << "\n                              computed    "<< H.determinant();
  100.  
  101.   cout << "\nDone\n";
  102. }
  103.  
  104. /*
  105.  *------------------------------------------------------------------------
  106.  *            Verify matrix multiplications
  107.  */
  108.  
  109. static void test_mm_multiplications(const int msize)
  110. {
  111.   cout << "\n---> Verify matrix multiplications\n"
  112.           "for matrices of the characteristic size " << msize << "\n\n";
  113.  
  114.   {
  115.     cout << "Test inline multiplications of the UnitMatrix" << endl;
  116.     Matrix m(-1,msize,-1,msize);
  117.     Matrix u(Matrix::Unit,m);
  118.     m.hilbert_matrix(); m(3,1) = M_PI;
  119.     u *= m;
  120.     verify_matrix_identity(u,m);
  121.   }
  122.  
  123.   {
  124.     cout << "Test inline multiplications by a DiagMatrix" << endl;
  125.     Matrix m(msize+3,msize);
  126.     m.hilbert_matrix(); m(1,3) = M_PI;
  127.     Vector v(msize);
  128.     register int i;
  129.     for(i=v.q_lwb(); i<=v.q_upb(); i++)
  130.       v(i) = 1+i;
  131.     Matrix diag(msize,msize);
  132.     (MatrixDiag)diag = v;
  133.     Matrix eth = m;
  134.     for(i=eth.q_row_lwb(); i<=eth.q_row_upb(); i++)
  135.       for(register int j=eth.q_col_lwb(); j<=eth.q_col_upb(); j++)
  136.     eth(i,j) *= v(j);
  137.     m *= diag;
  138.     verify_matrix_identity(m,eth);
  139.   }
  140.  
  141.   {
  142.     cout << "Test XPP = X where P is a permutation matrix\n";
  143.     Matrix m(msize-1,msize);
  144.     m.hilbert_matrix(); m(2,3) = M_PI;
  145.     Matrix eth = m;
  146.     Matrix p(msize,msize);
  147.     for(register int i=p.q_row_lwb(); i<=p.q_row_upb(); i++)
  148.       p(p.q_row_upb()+p.q_row_lwb()-i,i) = 1;
  149.     m *= p;
  150.     m *= p;
  151.     verify_matrix_identity(m,eth);
  152.   }
  153.  
  154.   {
  155.     cout << "Test general matrix multiplication through inline mult" << endl;
  156.     Matrix m(msize-2,msize);
  157.     m.hilbert_matrix(); m(3,3) = M_PI;
  158.     Matrix mt(Matrix::Transposed,m);
  159.     Matrix p(msize,msize);
  160.     p.hilbert_matrix();
  161.     MatrixDiag(p) += 1;
  162.     Matrix mp(m,Matrix::Mult,p);
  163.     Matrix m1 = m;
  164.     m *= p;
  165.     verify_matrix_identity(m,mp);
  166.     Matrix mp1(mt,Matrix::TransposeMult,p);
  167.     verify_matrix_identity(m,mp1);
  168.     assert( !(m1 == m) );
  169.     Matrix mp2(Matrix::Zero,m1);
  170.     assert( mp2 == 0 );
  171.     mp2.mult(m1,p);
  172.     verify_matrix_identity(m,mp2);
  173.   }
  174.  
  175.   {
  176.     cout << "Check to see UU' = U'U = E when U is the Haar matrix" << endl;
  177.     const int order = 5;
  178.     const int no_sub_cols = (1<<order) - 5;
  179.     Matrix haar_sub = haar_matrix(5,no_sub_cols);
  180.     Matrix haar_sub_t(Matrix::Transposed,haar_sub);
  181.     Matrix hsths(haar_sub_t,Matrix::Mult,haar_sub);
  182.     Matrix hsths1(Matrix::Zero,hsths); hsths1.mult(haar_sub_t,haar_sub);
  183.     Matrix hsths_eth(Matrix::Unit,hsths);
  184.     assert( hsths.q_nrows() == no_sub_cols && hsths.q_ncols() == no_sub_cols );
  185.     verify_matrix_identity(hsths,hsths_eth);
  186.     verify_matrix_identity(hsths1,hsths_eth);
  187.     
  188.     Matrix haar = haar_matrix(5);
  189.     Matrix unit(Matrix::Unit,haar);
  190.     Matrix haar_t(Matrix::Transposed,haar);
  191.     Matrix hth(haar,Matrix::TransposeMult,haar);
  192.     Matrix hht(haar,Matrix::Mult,haar_t);
  193.     Matrix hht1 = haar; hht1 *= haar_t;
  194.     Matrix hht2(Matrix::Zero,haar); hht2.mult(haar,haar_t);
  195.     verify_matrix_identity(unit,hth);
  196.     verify_matrix_identity(unit,hht);
  197.     verify_matrix_identity(unit,hht1);
  198.     verify_matrix_identity(unit,hht2);
  199.   }
  200.   cout << "\nDone\n";
  201. }
  202.  
  203. static void test_vm_multiplications(const int msize)
  204. {
  205.   cout << "\n---> Verify vector-matrix multiplications\n"
  206.           "for matrices of the characteristic size " << msize << "\n\n";
  207.   {
  208.     cout << "Check shrinking a vector by multiplying by a non-sq unit matrix"
  209.          << endl;
  210.     Vector vb(-2,msize);
  211.     for(register int i=vb.q_lwb(); i<=vb.q_upb(); i++)
  212.       vb(i) = M_PI - i;
  213.     assert( vb != 0 );
  214.     Matrix mc(1,msize-2,-2,msize);    // contracting matrix
  215.     mc.unit_matrix();
  216.     Vector v1 = vb;
  217.     Vector v2 = vb;
  218.     v1 *= mc;
  219.     v2.resize_to(1,msize-2);
  220.     verify_matrix_identity(v1,v2);
  221.   }
  222.  
  223.   {
  224.     cout << "Check expanding a vector by multiplying by a non-sq unit matrix"
  225.          << endl;
  226.     Vector vb(msize);
  227.     for(register int i=vb.q_lwb(); i<=vb.q_upb(); i++)
  228.       vb(i) = M_PI + i;
  229.     assert( vb != 0 );
  230.     Matrix me(2,msize+5,1,msize);    // expanding matrix
  231.     me.unit_matrix();
  232.     Vector v1 = vb;
  233.     Vector v2 = vb;
  234.     v1 *= me;
  235.     v2.resize_to(v1);
  236.     verify_matrix_identity(v1,v2);
  237.   }
  238.  
  239.   {
  240.     cout << "Check general matrix-vector multiplication" << endl;
  241.     Vector vb(msize);
  242.     for(register int i=vb.q_lwb(); i<=vb.q_upb(); i++)
  243.       vb(i) = M_PI + i;
  244.     Matrix vm(msize,1);
  245.     MatrixColumn(vm,1) = vb;
  246.     Matrix m(0,msize,1,msize);
  247.     m.hilbert_matrix();
  248.     vb *= m;
  249.     assert( vb.q_lwb() == 0 );
  250.     Matrix mvm(m,Matrix::Mult,vm);
  251.     verify_matrix_identity(vb,mvm);
  252.   }
  253.  
  254.   cout << "\nDone\n";
  255. }
  256.  
  257. static void test_inversion(const int msize)
  258. {
  259.   cout << "\n---> Verify matrix inversion for square matrices\n"
  260.           "of size " << msize << "\n\n";
  261.   {
  262.     cout << "Test invesion of a diagonal matrix" << endl;
  263.     Matrix m(-1,msize,-1,msize);
  264.     Matrix mi(Matrix::Zero,m);
  265.     for(register int i=m.q_row_lwb(); i<=m.q_row_upb(); i++)
  266.       mi(i,i) = 1/(m(i,i) = i-m.q_row_lwb()+1);
  267.     Matrix mi1(Matrix::Inverted,m);
  268.     m.invert();
  269.     verify_matrix_identity(m,mi);
  270.     verify_matrix_identity(mi1,mi);
  271.   }
  272.  
  273.   {
  274.     cout << "Test invesion of an orthonormal (Haar) matrix" << endl;
  275.     Matrix m = haar_matrix(3);
  276.     Matrix morig = m;
  277.     Matrix mt(Matrix::Transposed,m);
  278.     double det = -1;        // init to a wrong val to see if it's changed
  279.     m.invert(&det);
  280.     assert( (det-1) <= FLT_EPSILON );
  281.     verify_matrix_identity(m,mt);
  282.     Matrix mti(Matrix::Inverted,mt);
  283.     verify_matrix_identity(mti,morig);
  284.   }
  285.  
  286.   {
  287.     cout << "Test invesion of a good matrix with diagonal dominance" << endl;
  288.     Matrix m(msize,msize);
  289.     m.hilbert_matrix();
  290.     MatrixDiag(m) += 1;
  291.     Matrix morig = m;
  292.     double det_inv = 0;
  293.     const double det_comp = m.determinant();
  294.     m.invert(&det_inv);
  295.     cout << "\tcomputed determinant             " << det_comp << endl;
  296.     cout << "\tdeterminant returned by invert() " << det_inv << endl;
  297.  
  298.     cout << "\tcheck to see M^(-1) * M is E" << endl;
  299.     Matrix mim(m,Matrix::Mult,morig);
  300.     Matrix unit(Matrix::Unit,m);
  301.     verify_matrix_identity(mim,unit);
  302.  
  303.     cout << "\tcheck to see M * M^(-1) is E" << endl;
  304.     Matrix mmi = morig; mmi *= m;
  305.     verify_matrix_identity(mmi,unit);
  306.   }
  307.     
  308.   cout << "\nDone\n";
  309. }
  310.  
  311.  
  312. /*
  313.  *------------------------------------------------------------------------
  314.  *                Root module
  315.  */
  316.  
  317. main()
  318. {
  319.   cout << "\n\n" << _Minuses << 
  320.           "\n\t\tVerify Advanced Operations on Matrices\n";
  321.  
  322.   test_determinant(20);
  323.   test_mm_multiplications(20);
  324.   test_vm_multiplications(20);
  325.  
  326.   test_inversion(20);
  327. }
  328.  
  329.